TODO
options(timeout = 10000)
dir.create("data/steinbock/raw", recursive = TRUE)
## Warning in dir.create("data/steinbock/raw", recursive = TRUE): 'data/steinbock/
## raw' already exists
download.file("https://zenodo.org/record/6642699/files/panel.csv",
"data/steinbock/panel.csv")
download.file("https://zenodo.org/record/5949116/files/Patient1.zip",
"data/steinbock/raw/Patient1.zip")
download.file("https://zenodo.org/record/5949116/files/Patient2.zip",
"data/steinbock/raw/Patient2.zip")
download.file("https://zenodo.org/record/5949116/files/Patient3.zip",
"data/steinbock/raw/Patient3.zip")
download.file("https://zenodo.org/record/5949116/files/Patient4.zip",
"data/steinbock/raw/Patient4.zip")
download.file("https://zenodo.org/record/5949116/files/compensation.zip",
"data/compensation.zip")
unzip("data/compensation.zip", exdir="data", overwrite=TRUE)
unlink("data/compensation.zip")
download.file("https://zenodo.org/record/5949116/files/sample_metadata.xlsx",
destfile = "data/sample_metadata.xlsx")
download.file("https://zenodo.org/record/7079294/files/gated_cells.zip",
"data/gated_cells.zip")
unzip("data/gated_cells.zip", exdir="data", overwrite=TRUE)
unlink("data/gated_cells.zip")
if (!requireNamespace("devtools", quietly = TRUE)) install.packages("devtools")
devtools::install_github(c("BodenmillerGroup/imcRtools",
"BodenmillerGroup/cytomapper"))
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("pheatmap", "viridis",
"tiff", "distill", "openxlsx", "ggrepel", "patchwork",
"mclust", "RColorBrewer", "uwot", "Rtsne", "caret",
"randomForest", "ggridges", "gridGraphics", "scales",
"CATALYST", "scuttle", "scater", "dittoSeq",
"tidyverse", "batchelor", "bluster","scran"))
steinbock preprocess imc images --hpf 50
## bash: no job control in this shell
## 2022-09-27 14:07:32,574 INFO steinbock - img/Patient4_005.tiff
## 2022-09-27 14:07:35,642 INFO steinbock - img/Patient4_006.tiff
## 2022-09-27 14:07:37,683 INFO steinbock - img/Patient4_007.tiff
## 2022-09-27 14:07:40,122 INFO steinbock - img/Patient4_008.tiff
## 2022-09-27 14:07:42,320 INFO steinbock - img/Patient3_001.tiff
## 2022-09-27 14:07:44,818 INFO steinbock - img/Patient3_002.tiff
## 2022-09-27 14:07:51,711 INFO steinbock - img/Patient3_003.tiff
## 2022-09-27 14:08:04,518 INFO steinbock - img/Patient2_001.tiff
## 2022-09-27 14:08:06,077 INFO steinbock - img/Patient2_002.tiff
## 2022-09-27 14:08:09,393 INFO steinbock - img/Patient2_003.tiff
## 2022-09-27 14:08:12,661 INFO steinbock - img/Patient2_004.tiff
## 2022-09-27 14:08:15,646 INFO steinbock - img/Patient1_001.tiff
## 2022-09-27 14:08:19,383 INFO steinbock - img/Patient1_002.tiff
## 2022-09-27 14:08:21,434 INFO steinbock - img/Patient1_003.tiff
## 2022-09-27 14:08:35,616 INFO steinbock - images.csv
The step took 1.96 minutes.
steinbock segment deepcell --minmax
## bash: no job control in this shell
## /opt/venv/lib/python3.8/site-packages/deepcell_toolbox/deep_watershed.py:179: FutureWarning: `selem` is a deprecated argument name for `h_maxima`. It will be removed in version 1.0. Please use `footprint` instead.
## markers = h_maxima(image=maxima,
## 2022-09-27 14:09:31,316 INFO steinbock - masks/Patient1_001.tiff
## 2022-09-27 14:09:50,100 INFO steinbock - masks/Patient1_002.tiff
## 2022-09-27 14:10:16,780 INFO steinbock - masks/Patient1_003.tiff
## 2022-09-27 14:10:35,263 INFO steinbock - masks/Patient2_001.tiff
## 2022-09-27 14:10:58,705 INFO steinbock - masks/Patient2_002.tiff
## 2022-09-27 14:11:26,376 INFO steinbock - masks/Patient2_003.tiff
## 2022-09-27 14:11:52,986 INFO steinbock - masks/Patient2_004.tiff
## 2022-09-27 14:12:23,620 INFO steinbock - masks/Patient3_001.tiff
## 2022-09-27 14:12:50,053 INFO steinbock - masks/Patient3_002.tiff
## 2022-09-27 14:13:19,963 INFO steinbock - masks/Patient3_003.tiff
## 2022-09-27 14:13:45,334 INFO steinbock - masks/Patient4_005.tiff
## 2022-09-27 14:14:10,888 INFO steinbock - masks/Patient4_006.tiff
## 2022-09-27 14:14:33,406 INFO steinbock - masks/Patient4_007.tiff
## 2022-09-27 14:14:57,807 INFO steinbock - masks/Patient4_008.tiff
The step took 6.48 minutes.
steinbock measure intensities
## bash: no job control in this shell
## 2022-09-27 14:15:13,624 INFO steinbock - intensities/Patient1_001.csv
## 2022-09-27 14:15:15,273 INFO steinbock - intensities/Patient1_002.csv
## 2022-09-27 14:15:18,603 INFO steinbock - intensities/Patient1_003.csv
## 2022-09-27 14:15:19,997 INFO steinbock - intensities/Patient2_001.csv
## 2022-09-27 14:15:21,412 INFO steinbock - intensities/Patient2_002.csv
## 2022-09-27 14:15:22,845 INFO steinbock - intensities/Patient2_003.csv
## 2022-09-27 14:15:26,029 INFO steinbock - intensities/Patient2_004.csv
## 2022-09-27 14:15:28,630 INFO steinbock - intensities/Patient3_001.csv
## 2022-09-27 14:15:30,687 INFO steinbock - intensities/Patient3_002.csv
## 2022-09-27 14:15:32,342 INFO steinbock - intensities/Patient3_003.csv
## 2022-09-27 14:15:33,638 INFO steinbock - intensities/Patient4_005.csv
## 2022-09-27 14:15:35,263 INFO steinbock - intensities/Patient4_006.csv
## 2022-09-27 14:15:37,000 INFO steinbock - intensities/Patient4_007.csv
## 2022-09-27 14:15:38,335 INFO steinbock - intensities/Patient4_008.csv
The step took 0.54 minutes.
steinbock measure regionprops
## bash: no job control in this shell
## 2022-09-27 14:15:44,211 INFO steinbock - regionprops/Patient1_001.csv
## 2022-09-27 14:15:46,308 INFO steinbock - regionprops/Patient1_002.csv
## 2022-09-27 14:15:48,631 INFO steinbock - regionprops/Patient1_003.csv
## 2022-09-27 14:15:50,725 INFO steinbock - regionprops/Patient2_001.csv
## 2022-09-27 14:15:52,862 INFO steinbock - regionprops/Patient2_002.csv
## 2022-09-27 14:15:54,442 INFO steinbock - regionprops/Patient2_003.csv
## 2022-09-27 14:15:58,003 INFO steinbock - regionprops/Patient2_004.csv
## 2022-09-27 14:16:00,138 INFO steinbock - regionprops/Patient3_001.csv
## 2022-09-27 14:16:02,029 INFO steinbock - regionprops/Patient3_002.csv
## 2022-09-27 14:16:03,929 INFO steinbock - regionprops/Patient3_003.csv
## 2022-09-27 14:16:05,437 INFO steinbock - regionprops/Patient4_005.csv
## 2022-09-27 14:16:07,950 INFO steinbock - regionprops/Patient4_006.csv
## 2022-09-27 14:16:09,749 INFO steinbock - regionprops/Patient4_007.csv
## 2022-09-27 14:16:11,344 INFO steinbock - regionprops/Patient4_008.csv
The step took 0.55 minutes.
steinbock measure neighbors --type expansion --dmax 4
## bash: no job control in this shell
## 2022-09-27 14:16:18,146 INFO steinbock - neighbors/Patient1_001.csv
## 2022-09-27 14:16:22,073 INFO steinbock - neighbors/Patient1_002.csv
## 2022-09-27 14:16:26,821 INFO steinbock - neighbors/Patient1_003.csv
## 2022-09-27 14:16:30,658 INFO steinbock - neighbors/Patient2_001.csv
## 2022-09-27 14:16:34,024 INFO steinbock - neighbors/Patient2_002.csv
## 2022-09-27 14:16:37,350 INFO steinbock - neighbors/Patient2_003.csv
## 2022-09-27 14:16:42,074 INFO steinbock - neighbors/Patient2_004.csv
## 2022-09-27 14:16:47,165 INFO steinbock - neighbors/Patient3_001.csv
## 2022-09-27 14:16:51,557 INFO steinbock - neighbors/Patient3_002.csv
## 2022-09-27 14:16:56,411 INFO steinbock - neighbors/Patient3_003.csv
## 2022-09-27 14:16:59,925 INFO steinbock - neighbors/Patient4_005.csv
## 2022-09-27 14:17:05,593 INFO steinbock - neighbors/Patient4_006.csv
## 2022-09-27 14:17:10,029 INFO steinbock - neighbors/Patient4_007.csv
## 2022-09-27 14:17:13,830 INFO steinbock - neighbors/Patient4_008.csv
The step took 1.04 minutes.
library(imcRtools)
spe <- read_steinbock("data/steinbock/")
The step took 0.92 minutes.
library(openxlsx)
library(tidyverse)
colnames(spe) <- paste0(spe$sample_id, "_", spe$ObjectNumber)
meta <- read.xlsx("data/sample_metadata.xlsx")
spe$patient_id <- as.vector(str_extract_all(spe$sample_id, "Patient[1-4]",
simplify = TRUE))
spe$ROI <- as.vector(str_extract_all(spe$sample_id, "00[1-8]",
simplify = TRUE))
spe$indication <- meta$Indication[match(spe$patient_id, meta$Sample.ID)]
rowData(spe)$use_channel <- !grepl("DNA|Histone", rownames(spe))
assay(spe, "exprs") <- asinh(counts(spe)/1)
The step took 0.03 minutes.
library(cytomapper)
images <- loadImages("data/steinbock/img/")
channelNames(images) <- rownames(spe)
The step took 0.27 minutes.
masks <- loadImages("data/steinbock/masks/", as.is = TRUE)
## All files in the provided location will be read in.
The step took 0 minutes.
patient_id <- str_extract_all(names(images), "Patient[1-4]", simplify = TRUE)
indication <- meta$Indication[match(patient_id, meta$Sample.ID)]
mcols(images) <- mcols(masks) <- DataFrame(sample_id = names(images),
patient_id = patient_id,
indication = indication)
The step took 0 minutes.
sce <- readSCEfromTXT("data/compensation/")
## Spotted channels: Y89, In113, In115, Pr141, Nd142, Nd143, Nd144, Nd145, Nd146, Sm147, Nd148, Sm149, Nd150, Eu151, Sm152, Eu153, Sm154, Gd155, Gd156, Gd158, Tb159, Gd160, Dy161, Dy162, Dy163, Dy164, Ho165, Er166, Er167, Er168, Tm169, Er170, Yb171, Yb172, Yb173, Yb174, Lu175, Yb176
## Acquired channels: Ar80, Y89, In113, In115, Xe131, Xe134, Ba136, La138, Pr141, Nd142, Nd143, Nd144, Nd145, Nd146, Sm147, Nd148, Sm149, Nd150, Eu151, Sm152, Eu153, Sm154, Gd155, Gd156, Gd158, Tb159, Gd160, Dy161, Dy162, Dy163, Dy164, Ho165, Er166, Er167, Er168, Tm169, Er170, Yb171, Yb172, Yb173, Yb174, Lu175, Yb176, Ir191, Ir193, Pt196, Pb206
## Channels spotted but not acquired:
## Channels acquired but not spotted: Ar80, Xe131, Xe134, Ba136, La138, Ir191, Ir193, Pt196, Pb206
assay(sce, "exprs") <- asinh(counts(sce)/5)
The step took 0.09 minutes.
plotSpotHeatmap(sce)
sce2 <- binAcrossPixels(sce, bin_size = 10)
The step took 0.37 minutes.
library(CATALYST)
bc_key <- as.numeric(unique(sce$sample_mass))
bc_key <- bc_key[order(bc_key)]
sce <- assignPrelim(sce, bc_key = bc_key)
## Debarcoding data...
## o ordering
## o classifying events
## Normalizing...
## Computing deltas...
sce <- estCutoffs(sce)
sce <- applyCutoffs(sce)
sce <- filterPixels(sce, minevents = 40, correct_pixels = TRUE)
The step took 0.81 minutes.
sce <- computeSpillmat(sce)
sm <- metadata(sce)$spillover_matrix
The step took 0.07 minutes.
library(dittoSeq)
library(patchwork)
rowData(spe)$channel_name <- paste0(rowData(spe)$channel, "Di")
isotope_list <- CATALYST::isotope_list
isotope_list$Ar <- 80
spe <- compCytof(spe, sm,
transform = TRUE, cofactor = 1,
isotope_list = isotope_list,
overwrite = FALSE)
before <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
assay.x = "exprs", assay.y = "exprs") +
ggtitle("Before compensation")
after <- dittoScatterPlot(spe, x.var = "Ecad", y.var = "CD303",
assay.x = "compexprs", assay.y = "compexprs") +
ggtitle("After compensation")
before + after
assay(spe, "counts") <- assay(spe, "compcounts")
assay(spe, "exprs") <- assay(spe, "compexprs")
assay(spe, "compcounts") <- assay(spe, "compexprs") <- NULL
channelNames(images) <- rowData(spe)$channel_name
adapted_sm <- adaptSpillmat(sm, paste0(rowData(spe)$channel, "Di"),
isotope_list = isotope_list)
## Compensation is likely to be inaccurate.
## Spill values for the following interactions
## have not been estimated:
## Ir191Di -> Ir193Di
## Ir193Di -> Ir191Di
images_comp <- compImage(images, adapted_sm)
plotPixels(images[5], colour_by = "Yb173Di",
image_title = list(text = "Yb173 (Ecad) - before",
position = "topleft"),
legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
plotPixels(images[5], colour_by = "Yb174Di",
image_title = list(text = "Yb174 (CD303) - before",
position = "topleft"),
legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
plotPixels(images_comp[5], colour_by = "Yb173Di",
image_title = list(text = "Yb173 (Ecad) - after",
position = "topleft"),
legend = NULL, bcg = list(Yb173Di = c(0, 4, 1)))
plotPixels(images_comp[5], colour_by = "Yb174Di",
image_title = list(text = "Yb174 (CD303) - after",
position = "topleft"),
legend = NULL, bcg = list(Yb174Di = c(0, 4, 1)))
channelNames(images_comp) <- rownames(spe)
The step took 12.53 minutes.
set.seed(20220118)
img_ids <- sample(seq_len(length(images_comp)), 3)
cur_images <- images_comp[img_ids]
cur_images <- normalize(cur_images, separateImages = TRUE)
cur_images <- normalize(cur_images, inputRange = c(0, 0.2))
plotPixels(cur_images,
mask = masks[img_ids],
img_id = "sample_id",
missing_colour = "white",
colour_by = c("CD163", "CD20", "CD3", "Ecad", "DNA1"),
colour = list(CD163 = c("black", "yellow"),
CD20 = c("black", "red"),
CD3 = c("black", "green"),
Ecad = c("black", "cyan"),
DNA1 = c("black", "blue")),
image_title = NULL,
legend = list(colour_by.title.cex = 0.9,
colour_by.labels.cex = 0.9))